#include "definitions.h"
#include "base_uday_wrapper.hpp"

#include "tbb/tick_count.h"

extern "C" {
  void heat_3d_uday_omp(double (*M)[2][Npadded][Npadded][Npadded], int pT, int pN);
  void heat_3d_uday_no_skew(double (*M)[2][Npadded][Npadded][Npadded], int pT, int pN);
  void heat_3d(double (*M)[2][Npadded][Npadded][Npadded], int pT, int pN);

  void heat_3d(double (*M)[2][Npadded][Npadded][Npadded], int pT, int pN) {
    int i, j, k, t;
    for (t=0; t<pT-1; t++) {
      for (i=1; i<pN+1; i++) {
        for (j=1; j<pN+1; j++) {
          for (k=1; k<pN+1; k++) {
            M[0][(t+1)%2][i][j][k] = 
              0.125 * (M[0][t%2][i+1][j][k] - 2.0 * M[0][t%2][i][j][k] + M[0][t%2][i-1][j][k])
              + 0.125 * (M[0][t%2][i][j+1][k] - 2.0 * M[0][t%2][i][j][k] + M[0][t%2][i][j-1][k])
              + 0.125 * (M[0][t%2][i][j][k-1] - 2.0 * M[0][t%2][i][j][k] + M[0][t%2][i][j][k+1])
              + M[0][t%2][i][j][k];
          }
        }
      }
    }
  }
}

void Uday::operator()() {
   tbb::tick_count tick = tbb::tick_count::now();

   if (version == 1)
     heat_3d_uday_omp((double (*)[2][Npadded][Npadded][Npadded]) M, pT, pN);  
   else {
     assert(version == 2);
     heat_3d_uday_no_skew((double (*)[2][Npadded][Npadded][Npadded]) M, pT, pN);  
   }

   printf("Total time for init and parallel execution of %ld iterations of a" \
          "stencil of size %ld x %ld x %ld: %g sec\n", pT, pN, pN, pN, (tbb::tick_count::now()-tick).seconds());
}

